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ABSTRACT 

We investigate the signatures left by massive neutrinos on the spatial distribution of neutral hy¬ 
drogen (HI) in the post-reionization era by running hydrodynamic simulations that include massive 
neutrinos as additional collisionless particles. We find that halos in massive/massless neutrino cos¬ 
mologies host a similar amount of neutral hydrogen, although for a fixed halo mass, on average, the 
HI mass increases with the sum of the neutrino masses. Our results show that HI is more strongly 
clustered in cosmologies with massive neutrinos, while its abundance, Hhi(^), is lower. These effects 
arise mainly from the impact of massive neutrinos on cosmology: they suppress both the amplitude 
of the matter power spectrum on small scales and the abundance of dark matter halos. Modelling the 
HI distribution with hydrodynamic simulations at z > 3, and a simple analytic model at z < 3, we use 
the Fisher matrix formalism to conservatively forecast the constraints that Phase 1 of the Square Kilo¬ 
metre Array (SKA) will place on the sum of neutrino masses. Mi, = We find that with 10,000 

hours of interferometric observations at 3 z < 6 from a deep and narrow survey with SKAl-LOW, 
the sum of the neutrino masses can be measured with an error a{M„) < 0.3 eV (95% CL). Similar 
constraints can be obtained with a wide and deep SKAl-MID survey at z < 3, using the single-dish 
mode. By combining data from MID, LOW, and Planck, plus priors on cosmological parameters from 
a Stage IV spectroscopic galaxy survey, the sum of the neutrino masses can be determined with an 
error a{M„) ~ 0.06 eV (95% CL). 

Subject headings: massive neutrinos, intensity mapping, cosmology 


1. INTRODUCTION 

The standard model of particle physics describes neu¬ 
trinos as neutral spin-1/2, massless fermions, organized 
into three families with three different flavours: Ug, 

It has been observed that neutrinos can change their 
flavour as they propagate through space, however. This 
phenomenon, known as neutrino oscillations, implies that 
neutrinos are massive. Measurements of the neutrino os¬ 
cillations from laboratory experiments have allowed us to 
estimate the mass-square differences among the different 
neutrino mass eigenstates to be (Fogli et al. 2012; Forero 
et al. 2012): 

A mi2 = 7.5 X lO"'^ eV^ (1) 

I A 771^31= 2.3 X 10-3 eV^ (2) 

which implies that at least two of the three neutrino fam¬ 
ilies are massive. A lower bound on the sum of the neu¬ 
trino masses can be set from the above measurements: 
Ml, = ^ 0-06 eV. Unfortunately, the above con¬ 

straints do not allow us to determine which neutrino is 
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the lightest, or whether it is massless or massive. This 
gives rise to two different hierarchies: a normal hierar¬ 
chy in which 0 ^ toi < 7712 < 7713 , and an inverted hierarchy 
where 0 ^ 7773<7711 < 777 , 2 . 

The fact that neutrinos are massive is one of the clear¬ 
est indications of physics beyond the particle physics 
standard model, and so two of the most important ques¬ 
tions in modern physics are: (a) what are the masses of 
the neutrinos, and (b) which hierarchy do they conform 
to? Answering these questions with laboratory exper¬ 
iments is extremely challenging. For instance, current 
bounds on the mass of the electron anti-neutrino from the 
KATRIN^ experiment are m{Dg)<2.3 eV (Kraus et al. 
2005) and are expected to improve to m{ve)<0.2 eV in 
the coming years. 

On the other hand, tight upper limits on the neu¬ 
trino masses have already been obtained by using cos¬ 
mological observables such as the anisotropies in the 
cosmic microwave background (GMB), the clustering 
of galaxies, the abundance of galaxy clusters, the dis¬ 
tortion in the shape of galaxies by weak leasing, the 
Lya forest, etc. (Hannestad 2003; Reid et al. 2010; 
Thomas et al. 2010; Swanson et al. 2010; Saito et al. 

^ https://www.katriii.kit.edu/ 
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2011; Abazajian et al. 2011; de Putter et al. 2012; Xia 
et al. 2012; Rienier-S0rensen et al. 2012; Zhao et al. 2012; 
Costanzi Alunno Cerbolini et al. 2013; Basse et al. 2013; 
Planck Collaboration 2013; Costanzi et al. 2013, 2014; 
Wyman et al. 2014; Battye & Moss 2014; Hamann & 
Hasenkamp 2013; Beutler et al. 2014; Giusarma et al. 
2014; Palanque-Delabrouille et al. 2015b; Planck Collab¬ 
oration 2015; Palanque-Delabrouille et al. 2015a). These 
constraints arise because massive neutrinos delay the 
matter-radiation equality time and slow down the growth 
of matter perturbation on small scales. At linear order, 
the resulting effects on the CMB and matter power spec¬ 
trum are well known and understood, making cosmolog¬ 
ical observables extremely useful tools for putting upper 
limits on the sum of the neutrino masses. 

Currently, the tightest upper limit on the sum of the 
neutrino masses comes from combining data from the 
CMB, baryonic acoustic oscillations (BAO) and the Lya 
forest: M^, < 0.12 eV (95% CL) (Palanque-Delabrouille 
et al. 2015a). These limits also have strong implications 
for particle physics experiments like neutrinoless double 
beta decay (Dell’Oro et al. 2015). 

A new cosmological observable has recently been pro¬ 
posed that is expected to play an important role in future 
cosmology: 21cm intensity mapping (Bharadwaj et al. 
2001; Bharadwaj & Sethi 2001; Battye et al. 2004; Mc- 
Quinn et al. 2006; Chang et al. 2008; Loeb & Wyithe 
2008; Bull et al. 2015). The idea of this technique is to 
measure the integrated 21cm emission from unresolved 
galaxies by performing a low angular resolution survey 
(Santos et al. 2015). Since neutral hydrogen (HI) is a 
tracer of the underlying matter distribution of the Uni¬ 
verse on large scales, the HI power spectrum is expected 
to follow the shape of the matter power spectrum, but 
with a different amplitude (the HI bias). This should 
allow tight constraints to be placed on cosmological pa¬ 
rameters through measurements of the power spectrum 
of the 21cm field (Bull et al. 2015). 

Intensity mapping therefore constitutes a promis¬ 
ing new cosmological observable that can be used to 
constrain the neutrino masses (Loeb & Wyithe 2008; 
Pritchard & Pierpaoli 2008; Tegmark & Zaldarriaga 
2009; Metcalf 2010; Abazajian et al. 2011; Oyama et al. 
2013; Shimabukuro et al. 2014). In order to do that, 
however, one needs to understand how the 21cm power 
spectrum is affected by the presence of massive neutri¬ 
nos. The aim of this paper is to investigate the signatures 
left by massive neutrinos on the 21cm power spectrum 
in the post-reionization universe, in both the linear and 
fully non-linear regimes. 

We begin by studying the effects that massive neutri¬ 
nos have on the spatial distribution of neutral hydrogen 
in real-space. We do this by running hydrodynamic sim¬ 
ulations with massless and massive neutrinos. We inves¬ 
tigate how the presence of massive neutrinos affects the 
HI abundance and clustering properties and, ultimately, 
the signatures left by neutrinos in the 21cm power spec¬ 
trum. 

We also forecast the constraints that the future Square 
Kilometre Array (SKA) radio telescope will place on the 
sum of the neutrino masses. We do this using the Fisher 
matrix formalism, where the spatial distribution of neu¬ 
tral hydrogen is modeled using hydrodynamic simula¬ 
tions at redshifts 3 ^ ^ 5.5 (a redshift range where 


SKAl-LOW will collect data), and with a simple analytic 
model at redshifts z ^ 3 (the redshift range covered by 
SKAl-MID). Our approach is conservative in the sense 
that: (i) at z > 3 we compare models of the HI distri¬ 
bution using 4 different methods; and (ii) we embed the 
information from the 21cm power spectrum in the Fisher 
matrix forecasts in a conservative way. 

This paper is organized as follows. In Sec. 2 we de¬ 
scribe the set of hydrodynamic simulations carried out 
for this work. Our simulations do not account for two 
crucial processes needed to properly model the spatial 
distribution of neutral hydrogen: HI self-shielding, and 
the formation of molecular hydrogen. We correct the out¬ 
puts of our simulations a posteriori to account for these 
effects, depicting the four different methods we use to 
achieve this in Sec. 3. We also describe the method we 
use to model the spatial distribution of neutral hydrogen 
at redshifts z < 3, which is not covered by our hydro- 
dynamic simulations. In Sec. 4 we investigate the effect 
of massive neutrinos on the abundance and spatial dis¬ 
tribution of neutral hydrogen. We present our forecasts 
on the neutrino masses in Sec. 5 and, finally, draw the 
main conclusions of this paper in Sec. 6. 

2 . HYDRODYNAMIC SIMULATIONS 

We model the spatial distribution of neutral hydro¬ 
gen by running high-resolution hydrodynamic simula¬ 
tions in 14 different cosmological models. The values 
of the cosmological parameters of our fiducial model are 
= 0.3175, Hcdm = 0.2865, Ob = 0.049, = 0, 

Oa = 0.6825, h = 0.6711, n, = 0.9624 and CTg = 0.834, 
in excellent agreement with the latest results from Planck 
(Planck Collaboration 2015). In all simulations we have 
assumed a flat cosmology, and therefore the value of Oa 
is set to 1 — Om, with 0^ = Ocdm -I- Ob -I- with 
Oj^/i^ = My/(94.1 eV). A summary of our simulation 
suite is shown in Table 1. 

The simulations can be split into two different groups. 
On one hand we have simulations in which the value of 
one of the parameters, Ocdm, Ob, O^, h, rig, Ag, is varied 
(with respect to the value in the fiducial model) while 
the values of the other parameters are kept fixed. We 
use these simulations in our Fisher matrix analysis to in¬ 
vestigate degeneracies between cosmological parameters, 
and to forecast the constraints that the SKA will place on 
the neutrino masses. The simulations belonging to this 
group are (J",C+,C“,H+,"H”, A”*", A“). 

On the other hand we have simulations in which we 
vary the value of O,^, but keep O^ fixed. This is the most 
natural choice to investigate the effect of massive neutri¬ 
nos, since we assume that a fraction of the total matter 
content of the Universe is made up of neutrinos. We have 
run one simulation with = 0.3 eV, and another with 
Ml, = 0.6 eV. Even though these neutrino masses are 
ruled out by the most recent constraints that combine 
CMB, BAO and Lya-forest data, our purpose in this pa¬ 
per is to investigate the impact of neutrino masses on the 
HI spatial distribution. Note that for a realistic sum of 
neutrino masses, the effect will be very small and may be 
completely hidden by sample variance in our simulations. 
Therefore, we decided to run simulations with neutrino 
masses higher than current bounds to properly resolve 
the effects of neutrinos on the HI distribution. For the 
sum of the neutrino masses we use to run the Simula- 
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Name 

Box 

(/i~^Mpc) 

^cdm 

Db 

flu 

Da 


h 

ris 

10 ^ As 

^■8,0 

T 

50 

0.2685 

0.049 

0.0 

0.6825 

0 

0.67 

0.9624 

2.13 

0.834 

I/+ 

50 

0.2685 

0.049 

0.007075 

0.675425 

0 

0.67 

0.9624 

2.13 

0.778 


50 

0.261425 

0.049 

0.007075 

0.6825 

0 

0.67 

0.9624 

2.13 

0.764 

iy + + 

50 

0.25435 

0.049 

0.01415 

0.6825 

0 

0.67 

0.9624 

2.13 

0.693 

c+ 

50 

0.287 

0.049 

0.0 

0.664 

0 

0.67 

0.9624 

2.13 

0.868 

c- 

50 

0.25 

0.049 

0.0 

0.701 

0 

0.67 

0.9624 

2.13 

0.797 

B+ 

50 

0.2685 

0.055 

0.0 

0.6765 

0 

0.67 

0.9624 

2.13 

0.816 

B- 

50 

0.2685 

0.043 

0.0 

0.6885 

0 

0.67 

0.9624 

2.13 

0.853 

H+ 

50 

0.2685 

0.049 

0.0 

0.6825 

0 

0.71 

0.9624 

2.13 

0.886 

n- 

50 

0.2685 

0.049 

0.0 

0.6825 

0 

0.63 

0.9624 

2.13 

0.777 


50 

0.2685 

0.049 

0.0 

0.6825 

0 

0.67 

1.0009 

2.13 

0.846 

M- 

50 

0.2685 

0.049 

0.0 

0.6825 

0 

0.67 

0.9239 

2.13 

0.822 

A+ 

50 

0.2685 

0.049 

0.0 

0.6825 

0 

0.67 

0.9624 

2.45 

0.894 

A- 

50 

0.2685 

0.049 

0.0 

0.6825 

0 

0.67 

0.9624 

1.81 

0.769 


TABLE 1 

Summary of our simulation suite. The simulation name indicates the parameter that has been varied with respect to the 
FIDUCIAL, T, MODEL. ThE SUPERSCRIPT, +/“, DESIGNATES WHETHER THE VARIATION IS POSITIVE/nEGATIVE. ThE SIMULATIONS AND 
Um'*' are SIMULATIONS WITH MASSIVE NEUTRINOS HAVING A VALUE OF EQUAL TO THE ONE OF THE FIDUCIAL MODEL. 


tions, 0.3 eV and 0.6 eV, the neutrino masses are almost 
perfectly degenerate, so there is no need to distinguish 
between the three different families. We use these simu¬ 
lations to investigate the impact of massive neutrinos on 
the matter and HI spatial distribution. The simulations 
belonging to this group are (J^, 1 ^+, i/++). 

In each simulation we follow the evolution of 512^ CDM 
and 512^ baryon particles (plus 512^ neutrino particles 
for simulations with fli, > 0) in a periodic box of size 50 
comoving /i“^Mpc, down to redshift 3. For each simula¬ 
tion we save snapshots at redshifts 5.5, 5, 4.5, 4, 3.5 and 
3. Note that evolving our simulations down to z = 0, 
for all the cosmological models considered in this pa¬ 
per, would be extremely computationally expensive, so 
at redshifts lower than z = 3 we model the spatial distri¬ 
bution of neutral hydrogen using a simple analytic model 
described in Sec. 3.5. 

The simulations were run using the TreePM-|-SPH 
code GADGET-III (Springel 2005). They incorporate 
radiative cooling by hydrogen and helium, as well as 
heating by a uniform UV background. Both the cool¬ 
ing routine and the UV background have been modi¬ 
fied to obtain the desired thermal history, which cor¬ 
responds to the reference model of Viel et al. (2013) 
that has been shown to provide a good fit to the sta¬ 
tistical properties of the transmitted Lyman-a flux. In 
our simulations, hydrogen reionization takes place^ at 
z ~ 12 and the temperature-density relation for the low- 
density IGM T = To(z)(l -I- has 7 (z) = 1.3 and 

^ Note that this reionization redshift is in slight tension with 
the latest Planck results (Planck Collaboration 2015; Mitra et al. 
2015). We do not expect our conclusions to be affected by this 
however, since we use the same reionization history for all models, 
and are only interested in studying relative effects. 


To(z = 2.4,3,4) = (16500,15000,10000) K. For every 
simulation we generate 5000 quasar mock spectra at the 
snapshot redshifts, and tune the strength of the UV back¬ 
ground to reproduce the observed mean transmitted flux 
of the Lyman-a forest; this information is needed for 
two of the HI modeling methods (the pseudo-RT 1 and 
pseudo-RT 2 methods, described below). Star forma¬ 
tion is modeled using the multi-phase effective model of 
Springel & Hernquist (2003). 

The simulation initial conditions are generated at z = 
99 using the Zel’dovich approximation. We compute 
the transfer functions of the different components using 
CAMB (Lewis et al. 2000). In simulations with massive 
neutrinos, the initial conditions were generated taking 
into account the scale-dependent growth present in those 
cosmological models. We note that the random seeds 
used to generate the initial conditions are the same in all 
simulations. 

We identify dark matter halos using both the Friends- 
of-Friends (FoF) algorithm (Davis et al. 1985) with b = 
0.2 and SUBFIND (Springel et al. 2001; Dolag et al. 
2009). We require that a minimum of 32 CDM particles 
belong to the FoF halo to identify it. 


3. HI DISTRIBUTION 

Our hydrodynamic simulations do not take into ac¬ 
count two crucial physical processes needed to properly 
simulate the spatial distribution of neutral hydrogen: 
the formation of molecular hydrogen (H 2 ), and HI self¬ 
shielding. In this section we describe the various meth¬ 
ods we use to correct for those two processes. We con¬ 
sider four methods here: two models (pseudo-RT 1 and 
pseudo-RT 2) that aim to mimic the result of a full ra- 
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Pseudo-RT 1 

(fiducial method) 

Method 

Pseudo-RT 2 

Halo-based 1 

Halo-based 2 

HI self-shielding 

✓ 

✓ 

✓ 

✓ 

Molecular hydrogen 

✓ 

✓ 

X 

X 

Modeling 2 ) function 

X 

X 

✓ 

✓ 

HI assigned to all gas particles 

✓ 

✓ 

X 

X 

HI assigned only to gas particles in halos 

X 

X 

✓ 

✓ 

Value of r2Hi(-2) fixed a priori 

X 

X 

✓ 

✓ 

Reference work 

Rahmati et al. (2013a) 

Dave et al. (2013) 

Bagla et al. (2010) 

This work 


TABLE 2 

Differences between the four different methods used to model the spatial distribution oe neutral hydrogen. 


diative transfer calculation^, and two (halo-based 1 and 
halo-based 2) that were constructed by taking into ac¬ 
count the fact that all HI should be within dark matter 
halos (Villaescusa-Navarro et al. 2014b). These models 
have the objective of providing the shape and amplitude 
of the function Mhi(M, 2 ) (see subsection 3.5 for fur¬ 
ther details). We summarize the main features of each 
method in Table 2. 

In this section, we also explain the way we model the 
spatial distribution of neutral hydrogen at redshifts z<3. 
This redshift range is not covered by our simulations, but 
is needed to forecasts the constraints that the SKA will 
set on the neutrino masses. 


3.1. Pseudo-RT 1 (fiducial model) 

In the first pseudo-radiative transfer method, neutral 
hydrogen is assigned to every single particle in the sim¬ 
ulation. The HI self-shielding correction is modeled us¬ 
ing the fitting formula of Rahmati et al. (2013a) (see 
their appendix A) that was obtained by performing ra¬ 
diative transfer calculations on top of hydrodynamic sim¬ 
ulations. This correction states that the photo-ionization 
rate seen by a particular gas particle is a function of its 
density. 

The HI masses obtained in this way are further cor¬ 
rected to account for the presence of molecular hydro¬ 
gen, which is only assigned to star-forming particles. We 
assume that the ratio of molecular to neutral hydrogen 
scales with the pressure, P, as 


^H2 

Shi 



( 3 ) 


where Sh 2 and Shi are the molecular and neutral hy¬ 
drogen surface densities, respectively. In our analysis we 
use (Po,a) = (3.5 x 10^ cm“^K,0.5), where a is slightly 
different to the value measured by Blizt & Rosolowsky 
(Blitz & Rosolowsky 2006), a = 0.92. The reason for this 
choice is purely phenomenological - using a = 0.5, we ob- 


^ In both pseudo-RT methods, we only account for the radiation 
from the UV background and do not consider radiation from local 
sources (Miralda-Escude 2005; Schaye 2006; Rahmati et al. 2013b). 
We are interested only in studying relative differences here, rather 
than absolute quantities, and do not expect our conclusions to 
change by neglecting the radiation from local sources. 


tain much better agreement with the abundance of ab¬ 
sorbers with large column densities than with a = 0.92. 

We use the above procedure as our fiducial method 
to model the spatial distribution of neutral hydrogen in 
the various simulated cosmologies, and to investigate the 
effects induced by massive neutrinos on the spatial dis¬ 
tribution of neutral hydrogen (see Sec. 4). To study the 
robustness of our SKA neutrino masses forecasts, we also 
model the distribution of HI using three other methods, 
that we describe in the following subsections. 

3.2. Pseudo-RT 2 

The HI self-shielding correction can be implemented 
using a method proposed by Dave et al. (2013). The 
authors of that paper state that good agreement between 
this method and the full radiative transfer simulations 
by Faucher-Giguere et al. (2010) is achieved. It is also 
able to reproduce several observations, such as the HIMF 
(Haynes et al. 2011) sA, z = Q. We will briefly describe the 
method here, but refer the reader to Dave et al. (2013) 
for further details. 

For every gas particle in the simulation, the HI frac¬ 
tion in photo-ionization equilibrium is computed by tak¬ 
ing into account the strength of the UV background and 
the physical properties of the particle, such as its density 
and temperature. Then, the HI column density, for every 
gas particle, is computed by integrating the SPH kernel 
from the radius of the particle up to a given radius (if 
it exists), rthres, where it reaches a given threshold that 
we set to 10^^'^ cm“^. The method implements the HI 
self-shielded correction by assuming that 90% of the hy¬ 
drogen between r = 0 and r = rthres is fully neutral. 
Finally, the presence of molecular hydrogen is accounted 
for using the procedure described in the previous subsec¬ 
tion. 


3.3. Halo-based 1 

In contrast with the above two methods, where HI is 
assigned to each individual gas particle in the simulation 
according to its physical properties, the method depicted 
in this and in the next subsection are designed to assign 
HI to dark matter halos (see subsection 3.5 for further 
details). These methods assume that a dark matter halo 
of mass M at redshift z hosts an amount of HI given by 
the deterministic function Mhi(M, z). 
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Bagla et al. (2010) proposed a parametrization of the 
function Mhi(M, z) as 


Mhi(M, z) 


0 otherwise, 


( 4 ) 

where the values of the parameters M^i^{z) and Mi„ax(- 2 ) 
correspond to halos with circular velocities equal to 
Vmin = 30 km/s and Vmux = 200 km/s at redshift z, 
respectively. The value of fsiz) is tuned to reproduce 
the HI density parameter Hhi(z) (defined as the ratio 
between the comoving density of neutral hydrogen at red- 
shift z to the critical density at z = 0). For the redshifts 
covered by our simulations (3 ^ z ^ 5.5), we assume that 
the value of riHi(z:) does not depend on redshift, and that 
its value is equal to 10“^, both for halo-based 1 and for 
the halo-based 2 method. This is in excellent agreement 
with observations. 

In practice, this method works as follows: from a sim¬ 
ulation snapshot we identify all the dark matter halos. 
The total HI mass residing in a particular halo of mass 
M at redshift z is computed from Eq. 4. Finally, the 
HI mass within the halo is distributed according to some 
HI density profile, pm{M,z). We model the last step 
by splitting the total HI mass in a given halo equally 
amongst all the gas particles belonging to it. 

In Villaescusa-Navarro et al. (2014b), it was shown 
that the above method is capable of reproducing the 
damped Lyman alpha absorber (DLA) column density 
distribution function extremely well at redshifts z ^ 
[2 — 4]. In Padmanabhan et al. (2015) the authors showed 
that it also reproduces the HI bias at z = 0 and the prod¬ 
uct bm{z) X r2Hi(z) at z ~ 0.8 from 2Icm observations. 


3.4. Halo-based 2 

While the halo-based I model is capable of reproducing 
many observables, it does fail at reproducing one: the 
bias of the DLAs at z ~ 2.3 recently measured by the 
BOSS collaboration: 5dla = (2.17 ± 0.20)/3^-^^ (Font- 
Ribera et al. 2012), where Pf is the Lya forest redshift 
distortion parameter, whose value is of order I. Here we 
propose a simple model that is capable of reproducing 
this observable as well. We propose a functional form 


Mhi(M, z) 


fi{z)M if Mmi„(z) < M 
0 otherwise. 


where Mniin(z) is chosen to be the mass of dark matter 
halos with circular velocities equal to fniin = 62 km/s 
at redshift z, and / 4 (z) is a parameter whose value is 
set to reproduce the value of f2Hi(z). This simple model 
predicts a value of the HI bias at z = 2.3 (see Eq. 7) 
equal to 2.15, in perfect agreement with the observational 
measurements.'^ 


3.5. HIatz<S 

Running all of our high-resolution hydrodynamic sim¬ 
ulations down to z = 0 is infeasible given the computa¬ 
tional resources we have access to. On the other hand, 

^ Note that we are making the assumption that the bias of the 
DLAs is the bias of the HI. This assumption is reasonable since the 
amount of HI in Lyman limit systems and in the Lyo forest is of 
order ~ 10% 


the redshift range 0 ^ z<3 may be important when fore¬ 
casting the constraints on the neutrino masses that will 
be achievable with SKA. We therefore model the spatial 
distribution of neutral hydrogen at z<3 using a simple 
analytic model. 

In Villaescusa-Navarro et al. (2014b) it was shown that 
the fraction of HI outside dark matter halos in the post¬ 
reionization era is negligible. One can therefore make 
the assumption that all HI resides in dark matter halos. 
Under this assumption, and following the spirit of the 
halo model (Cooray & Sheth 2002), we can then predict 
the shape and amplitude of the HI power spectrum, in 
real-space, if we have the following ingredients: the halo 
mass function, n{M,z), the halo bias, b{M,z), the lin¬ 
ear matter power spectrum P/("(fc), and the functions 
z) and p-^pM, z). The functions Mhi(M, z) and 
PY{i{r\M,z) represent the average HI mass and density 
profile in a dark matter halo of mass M at redshift z. 

On large, linear, scales the HI power spectrum in real- 
space does not depend on the pHi(?’|Tf, z) function, but 
only on Mhi(M, z), and it is given by 

PBi{k,z) = bl^Pz)P^{k,z) , (6) 

where the HI bias, &hi(- 2) is given by 

7 , . , _ /o°°»^(^,2)&(^,^)^Hi(M,z)dM 

/“n(M,z)MHi(M,z)dM ’ ^ ^ 

The 2Icm power spectrum is given by 

P 2 ic^{k,z) = 5Tb{z)b^i{z) ^I -b ^/3(z) -b (8) 

X Pm (A:, z), 

where /3 is the redshift-space distortion parameter given 
by /3(z) = /(z)/6hi(z), with /(z) being the growth rate 
at redshift z. The third term on the right hand side 
arises from the Kaiser formula (Kaiser 1987). The value 
oi 5Ti,{z) is given by 

m{z) = 189 ^m[z)h mK , (9) 

where H{z) and Hq are the value of the Hubble param¬ 
eter at redshifts z and 0, respectively, h represents the 
value of Hq in units of 100 km s“^Mpc“^. Also, r2Hi(z) 
does not depend on piii{r\M, z), but only on the function 
Mhi(M,z): 

I r°° 

f2Hi(2:) =- / n{M,z)M}ii{M,z)dM . (10) 

Pc,0 Jo 

The above equations make clear the central role played 
by the function Myii{M,z) in studies related to 2Icm 
intensity mapping in the post-reionization era. Modeling 
this function is therefore all we need to predict the shape 
and amplitude of the HI/2Icm power spectrum on linear 
scales. 

We use the above equations to model the 2Icm power 
spectrum at redshifts z<3, with Mhi(M, z) given by the 
halo-based I prescription (Bagla et al. 2010). As dis¬ 
cussed in Sec. 3.3, the halo-based I model has three free 
parameters: Mmax{z) and / 3 (z). While the val¬ 

ues of the parameters Mmin{z) and Mmax{z) are chosen 
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to correspond to halos with circular velocities of 30 km/s 
and 200 km/s, the value of fsiz) is fixed by requiring 
that f2Hi(-2) reproduces the observational measurements. 
Since at z<3 observations disfavor models with constant 
r^Hij we follow Crighton et al. (2015) and assume a red- 
shift dependence nHi( 2 ) = 4x lO”"*)!which fully 
determines the value of fsiz). 

For the halo mass function and halo bias we use the 
Sheth & Tormen (1999) and Sheth, Mo & Tormen (Sheth 
et al. 2001) models, respectively. We use the matter 
power spectrum from halofit (Takahashi et al. 2012) 
to evaluate Pm{k,z) for a given cosmological model, 
and compute the Hl/21cm power spectra at redshifts 
z = {0,0.25,0.5,0.75,1,1.25,1.5,1.75,2,2.5}. 

In cosmologies with massive neutrinos we use the 
CDM-l-baryon field, instead of the total matter density 
field, to compute the 21cm power spectrum (last term on 
the r.h.s of Eq. 9), and to evaluate the halo mass func¬ 
tion and halo bias (Ichiki & Takada 2012; Castorina et al. 
2014). 

We now briefly discuss the differences between apply¬ 
ing the halo-based 1 model to our hydrodynamic simu¬ 
lations and the formalism we have described in this sub¬ 
section. By putting the HI in the gas particles belong¬ 
ing to the dark matter halos of our simulations, we not 
only model the function Mhi(M, z), but also the HI den¬ 
sity profile within halos, p}ii{r\M, z), which is ignored in 
the method above. The above formalism also implicitly 
assumes a scale-independent bias, and that the redshift- 
space distortions are accounted for by the Kaiser formula. 
The fully non-linear clustering and redshift-space distor¬ 
tions are taken into account by placing the HI in the sim¬ 
ulations, however. On large, linear scales, both methods 
should give the same results, while on small scales the 
method described in this section will break down. 

As such, we conclude by emphasizing that the above 
methodology is limited to linear scales, i.e. scales in 
which the HI bias is constant, and redshift-space distor¬ 
tions can be accounted for by using the Kaiser formula. 

4. EFFECT OF MASSIVE NEUTRINOS 

Here we study the effects induced by massive neutri¬ 
nos on the spatial distribution of neutral hydrogen. We 
investigate how the HI abundance and clustering proper¬ 
ties are affected by the presence of massive neutrinos by 
comparing the distribution of neutral hydrogen from the 
simulations {M^ = 0.3 eV) and {M^=0.6 eV) to 
the one from the fiducial simulation, JT (Mi^ = 0.0 eV). 

In Fig. 1 we show the spatial distribution of neutral 
hydrogen (top row), matter (middle row) and gas (bot¬ 
tom row) in a cosmology with massless neutrinos (left, 
simulation T) and with = 0.6 eV neutrinos (right, 
simulation The images have been created by tak¬ 

ing a slice of 2 h“^Mpc width. The spatial distribution 
of (total) matter is shown over the whole box (i.e. in a 
slice of 50 X 50 X 2 (h“^Mpc)^), while the gas and HI 
images display a zoom over the region marked with a 
red square. As can be seen, the differences in the spatial 
distribution of matter, gas, and (in particular) neutral 
hydrogen between the two models are very small. 

4.1. HI abundance 



Z 

Mo 

a 

(eV) 


(h-^MQ) 



3 

0.012 ±0.003 

0.699 ± 0.006 

0.0 

4 

0.015 ±0.003 

0.712 ±0.005 


5 

0.032 ± 0.008 

0.740 ± 0.006 


3 

0.015 ±0.005 

0.706 ±0.007 

0.3 

4 

0.031 ± 0.008 

0.732 ± 0.006 


5 

0.033 ± 0.006 

0.743 ± 0.005 


3 

0.016 ±0.005 

0.708 ±0.007 

0.6 

4 

0.027 ± 0.006 

0.728 ±0.006 


5 

0.061 ± 0.016 

0.760 ± 0.007 


TABLE 3 

Best-fit parameters for Mhi(M, z) = (M/Mo)“ as a 

FUNCTION OF REDSHIFT AND COSMOLOGY. 

We now investigate how the presence of massive neu¬ 
trinos impacts the function Mhi(M, z). By modeling the 
HI distribution using the pseudo-RT 1 method, we com¬ 
puted the HI mass within each dark matter halo in the 
simulations with massless and massive neutrinos. In the 
upper row of Fig. 2 we show the results at redshifts 
z = 3,4,5. 

For a hxed dark matter halo mass, we find that halos 
in the massless and massive neutrino models contain the 
same HI mass to well within one standard deviation, al¬ 
though halos in the massive neutrino cosmologies do tend 
to host a slightly higher HI mass. The HI mass excess 
with respect to the massless neutrino case is ~ 7% for 
the cosmology with 0.6 eV neutrinos, decreasing to ~ 3% 
for the 0.3 eV cosmology, with a very weak dependence 
on redshift. Our results show that the HI mass excess 
is not uniform in mass: the most massive halos host a 
higher fraction of HI compared with the low mass halos. 

The function Mhi (M, z) presents a cut-off at low 
masses, which can be seen clearly at z = 3. This cut-off 
is not physical, but is due to the resolution of our sim¬ 
ulations. In order to explore the physical cut-off arising 
from the fact that a minimum gas density and length is 
required to have self-shielded HI, simulations with higher 
resolution are needed. We leave this for a future work. 

For the mass range accessible in our simulations, we 
find that the Mhi ( A7, z) function can be fitted by a func¬ 
tion of the form: Mhi(M, z) = [M/Mo(z)]“^^^ In Table 
3 we show the best fit values of Mg and a for the three 
different cosmologies at redshifts z = 3,4,5. The value 
of the slope, a, increases with redshift for all of the mod¬ 
els, while at fixed redshift, a increases with the sum of 
the neutrino masses. This reflects a well known prop¬ 
erty: at a given redshift, the spatial distribution of mat¬ 
ter on small scales in a cosmology with massive neutrinos 
is effectively younger (has grown less) than its massless 
neutrino counterpart (Marulli et al. 2011; Villaescusa- 
Navarro et al. 2013a; Castorina et al. 2014; Costanzi 
et al. 2013; Villaescusa-Navarro et al. 2014a; Massara 
et al. 2014, 2015). 

Fig. 3 shows Hhi( 2 :) in each cosmology. We find that 
our fiducial model is capable of reproducing the values 
of Hhi(z) obtained from the observations of Songaila & 
Cowie (2010); Noterdaeme et al. (2012); Crighton et al. 
(2015) extremely well in the redshift range covered by 
our simulations. The redshift dependence of the func- 









Fig. 1. — Impact of massive neutrinos on the spatial distribution of neutral hydrogen (upper row), total matter (middle row) and gas 
(bottom row) at 2 : = 3. Panels on the left show the results for a massless neutrino cosmology while panels on the right are for a cosmological 
model with Mi, = 0.6 eV neutrinos. The middle panels display the spatial distribution of matter on in a slice of 50 X 50 X 2 (h“^Mpc)®, 
while top and bottom panels show a zoom into the region marked with a red square (the width of those slices is also 2 /i“^Mpc). 


tion r2Hi(.z) is very weak, although it is slightly more 
pronounced with increasing neutrino mass. 

For a fixed redshift, the value of decreases as 

the sum of the neutrino masses increases. This result 
can be understood if we look at Eq. 10 and take into 
account the fact that Mhi (M, z) barely changes among 
cosmologies with massive and massless neutrinos. The 
reason for the decrement in r2Hi(.z) is therefore due to 
the suppression in the abundance of halos that the pres¬ 
ence of massive neutrinos induces (Brandbyge et al. 2010; 
Marulli et al. 2011; Villaescusa-Navarro et al. 2013a; Cas- 
torina et al. 2014; Ichiki & Takada 2012; LoVerde 2014; 
Roncarelli et al. 2015; Castorina et al. 2015). Using our 


fitting function to Mhi(M, z) from Table 3 (with a cut¬ 
off at M = 2 x 10® we have checked that using 

the massive/massless neutrino halo mass functions repro¬ 
duces the decrement in Uhi(M, z) induced by neutrinos. 

We have also computed the HI column density distri¬ 
bution^ for the cosmologies with massless/massive neu¬ 
trinos; the results are shown on the bottom row of Fig. 2. 
At z = 3 we compare our results with the measure¬ 
ments by Noterdaeme et al. (2012). While these have 
a mean redshift (z) = 2.5, they cover a redshift range 

® See Appendix B of Villaescusa-Navarro et al. (2014b) for a 
description of the procedure used to calculate this. 











[ h -^ Mg] Mt,i„ [/i-i Mg] Mt^„ [/i-i Mg] 



Fig. 2.— Upper row: Function Mhi(M) at redshifts z = 3 (left), 2 = 4 (middle) and 2 = 5 (right) for the cosmological models with 
massless neutrinos (black), = 0.3 eV (magenta) and Mi, = 0.6 eV (green) when the HI is modeled using the pseudo-RT 1 method. For 
each dark matter halo we have computed the HI mass within it; the lines represent the running median, and error bars show the scatter 
around the mean. The bottom panels display the results normalized by the Mhi(M) function of the massless neutrino model. Bottom row: 
Same as above but for the HI column density distribution function. The observational measurements are from Noterdaeme et al. (2012) 
(2 = [2 — 3.5]) and Zafar et al. (2013) (using the whole redshift range: 2 = [1.5 — 5.0]) for the results at 2 = 3, and from Crighton et al. 
( 2015 ) (2 = [3.5 — 5.4]) and Zafar et al. (2013) (using only the redshift range 2 = [3.1 — 5.0]) for the plots at 2 = 4 and 2 = 5. 


z = [2 — 3.5], and we assume no redshift evolution down 
to z = 3. The abundance of DLAs in our simulations at 
z = 4 and z = 5 are compared against the recent mea¬ 
surements of Crighton et al. (2015), which have data in 
the redshift range z = [3.5 — 5.4]. The column density 
distribution function of sub-DLAs in our simulations is 
compared at z = 3 against the measurements by Zafar 
et al. (2013), obtained by using the whole redshift range 
z = [1.5 — 5.0], while at z = 4 and z = 5 we only use 
data in the redshift range z = [3.1 — 5.0]. We find that 
our fiducial model reproduces the observed abundance of 
DLAs and sub-DLAs very well at all redshifts. 

In the bottom panels we display the ratio between the 
HI column density distribution function of the models 
with massive and massless neutrinos. We find that mas¬ 
sive neutrinos suppress the abundance of DLAs and sub- 
DLAs at all redshifts, although the effect is stronger at 
higher redshift. This is due to the lower value of Dhi(2 :) 
that is present in cosmologies with massive neutrinos. 

We conclude that, for a given mass, dark matter ha¬ 


los in cosmologies with massless and massive neutrinos 
host, on average, the same amount of neutral hydrogen. 
The suppression on the halo mass function induced by 
massive neutrinos decreases the total amount of neutral 
hydrogen in the Universe, given the fact that only halos 
above a certain mass will host HI. This manifests itself in 
a lower value of Dhi(^:), and in a deficit in the abundance 
of DLAs and sub-DLAs in cosmologies with massive neu¬ 
trinos, with respect to the massless neutrinos model. 

4.2. HI clustering 

We now investigate the impact of massive neutrinos 
on the clustering properties of neutral hydrogen. We 
focus our attention in the HI power spectrum, Ppii{k,z), 
the HI bias, bpii{k,z), and the 2Icm power spectrum, 

.^21cm(^7 -^)- 

In the upper row of Fig. 4 we show the HI power spec¬ 
trum for the models with Mi, = 0.0,0.3, 0.6 eV neutrinos 
at redshifts z = 3,4, 5 when the HI is modeled using 
the pseudo-RT 1 method. Our results show that the HI 
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Fig. 3.— Redshift evolution of Qhi for cosmologies with mass¬ 
less neutrinos (black), = 0.3 eV (magenta) and = 0.6 eV 
(green) when the HI is modeled using the pseudo-RT 1 method. 
Values obtained from the observations of Noterdaeme et al. (2012); 
Crighton et al. (2015); Songaila & Cowie (2010) are shown with red, 
blue, and orange markers respectively. 


is more strongly clustered in cosmologies with massive 
neutrinos, with the clustering of the neutral hydrogen 
increasing with the sum of the neutrino masses.® The 
bottom panels of that row show the ratios of the HI power 
spectra for the massive and massless neutrino cosmolo¬ 
gies. We hnd that the increase in power in the massive 
neutrino cosmologies, relative to the fiducial model, is 
almost independent of scale. We find that differences in 
the HI clustering between models increase with redshift, 
however. At higher redshift, and for the Mi, = 0.6 eV 
model, the increase in power is also more scale-dependent 
than at lower redshift. We emphasize that the HI power 
spectrum, defined as Pm{k, z) = ((5Hi(fc, z)5^i{k, z)) does 
not depend on the value of Hhi(2 :), which, as we have seen 
above, is different for each cosmology. 

We can easily see why the HI is more strongly clus¬ 
tered in cosmologies with massive neutrinos if we take 
into account that, for a fixed dark matter halo mass, the 
halo bias increases with the sum of the neutrino masses 
(Marulli et al. 2011; Villaescusa-Navarro et al. 2014a; 
Castorina et al. 2014). This happens because, as we saw 
above, massive neutrinos induce a suppression of the halo 
mass function. Thus, for a given mass, halos in cosmolo¬ 
gies with massive neutrinos are rarer, and therefore they 
are more biased (for a detailed description see Massara 
et al. (2014)). Since Mhi(M, z) barely changes amongst 
the cosmologies (or, to be more precise, increases only 
slightly with the sum of the neutrino masses), we should 
therefore expect the neutral hydrogen to be more clus¬ 
tered in cosmologies with massive neutrinos, as we find.^ 

In the middle row of Fig. 4 we plot the HI bias, defined 
as hm{k,z) = PHi-m(fc, z)/T’m(A;, z), where Pm-m{k,z) 

® Note that neutral hydrogen is also more clustered in cosmolo¬ 
gies with warm dark matter than in the corresponding models with 
the standard cold dark matter, as shown in Carucci et al. (2015). 

^ We are also implicitly assuming that the low-mass cut-off in 
the Mhi(M, z) function is the same mass in cosmologies with mass¬ 
less/massive neutrinos, which we believe is a good assumption 
given the fact that that function barely changes amongst the cos¬ 
mologies in the mass range covered by our simulations. 


is the Hl-matter cross-power spectrum. As expected, the 
stronger clustering of HI in massive neutrino cosmologies 
is reflected in a higher HI bias in these cosmologies with 
respect to the fiducial, massless neutrinos model. On 
large scales, the HI bias increases by ^ 10% for the model 
with Ml, = 0.3 eV neutrinos, and ~ 30% for M„ = 0.6 
eV, with respect to the M^ = 0.0 eV cosmology. There 
is a very weak dependence on redshift. On large scales, 
the HI bias is more scale-dependent at higher redshift, 
as already noted by Carucci et al. (2015). 

It is also worth pointing out that at z = 3, while the 
HI bias is flat for wavenumbers k < 1 h Mpc~^ in the 
Ml, = 0 eV model, the models with massive neutrinos ex¬ 
hibit a HI bias that does depend mildly on scale. There 
are two explanations for this: 1) the bias is higher in 
cosmologies with massive neutrinos, and therefore more 
scale dependent; and 2) the HI bias, defined as above, 
is scale-dependent in massive neutrino cosmologies. The 
second option arises because it has recently been found 
that the halo bias in cosmologies with massive neutrinos 
is scale-dependent, even on very large scales (Villaescusa- 
Navarro et al. 2014a; Castorina et al. 2014). In Cas¬ 
torina et al. (2014) it was pointed out that if the halo 
bias is defined, in cosmologies with massive neutrinos, as 
P\i-ch(k)/Pcb{k), where cb denotes the cold dark mat¬ 
ter plus baryons field, then it becomes scale-independent 
on large scales and universal. We have checked whether 
defining the HI bias as 6 hi(^, z) = Pm-ch{k, z)/Pch{k, z) 
helps to decrease the scale-dependence of the bias at 
z = 3 for the cosmologies with massive neutrinos. The 
middle-left panel of Fig. 4 shows the results of defin¬ 
ing the HI bias with respect to the CDM-bbaryons field 
(dashed green line). We find that by using this deh- 
nition, the amplitude of the HI bias decreases because 
the CDM-l-baryon field is more strongly clustered that 
the total matter field in cosmologies with massive neu¬ 
trinos. We do not see a significant suppression of the 
scale-dependence, however. We therefore conclude that 
HI bias is more scale-dependent in cosmologies with mas¬ 
sive neutrinos because its value is higher (Scoccimarro 
et al. 2001; Cooray & Sheth 2002; Sefusatti & Scocci¬ 
marro 2005; Marfn et al. 2010). 

In the bottom row of Fig. 4 we show the 21cm power 
spectrum for the models with massless and massive neu¬ 
trinos. The amplitude of the 21cm power spectrum is 
lower in massive neutrino cosmologies than in the mass¬ 
less neutrinos model. While this result may seem surpris¬ 
ing, since we have seen above that HI clustering increases 
with the neutrino masses, it is straightforward to under¬ 
stand if we take into account that 

H21cm(fc) =^'(^)PHl(fc), (11) 

where P^i{k) denotes the HI power spectrum in redshift- 
space, and that STt{z) depends linearly on the value of 
r2Hi(z) (see Eq. 9). In other words, the amplitude of the 
21cm power spectrum not only depends on the HI clus¬ 
tering, but also on the value of Hhi(.z). Therefore, even 
if the HI is more clustered in cosmologies with massive 
neutrinos, the lower value of Hhi(z) in those cosmolo¬ 
gies drives the suppression of power we find in the 21cm 
power spectrum. 

In contrast with what happens with the HI power spec¬ 
trum, the ratio of the 21cm power spectra in the massive 
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Fig. 4.— Upper row: Power spectrum of the neutral hydrogen field for the fiducial model with massless neutrinos (black line), the model 
with Mu = 0.3 eV (magenta line) and the model with Mu = 0.6 eV (green) at 2 = 3 (left), 2 = 4 (middle), and 2 = 5 (right) when the HI 
is modeled using the pseudo-RT 1 method. The bottom panels display the ratio of the different HI power spectra to the fiducial model. 
The dashed vertical line represents the Nyquist frequency for our power spectrum measurements. Middle row: same as above but for the 
HI bias, computed as 5 hi(^) = ^Hi—m(^)/^m(^), with -PHi-m(^) being the Hl-matter cross-power spectrum. The green dashed line is HI 
bias computed with respect to the CDM-hbaryon field (see text for details). Bottom row: Same as above but for the 21cm field. 
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and massless neutrinos cosmologies exhibits a character¬ 
istic dependence on scale. On very small scales the ra¬ 
tio is flat, while it decreases around fc ~2 —3 /iMpc”^. 
The ratio appears to reach a minimum around k ^ 
0.3 /iMpc”^, although more realizations would be needed 
to confirm that this is not an artefact of sample variance. 

We conclude that the HI is more strongly clustered 
in cosmologies with massive neutrinos. This is because, 
for a fixed halo mass, the halo bias and HI mass increase 
with the sum of the neutrino masses. On the other hand, 
the amplitude of the 21cm power spectrum decreases as 
the neutrino masses increase because it depends on the 
total amount of neutral hydrogen, which is lower 

in massive neutrino cosmologies. 

5. SKA FORECASTS 

In this section we present Fisher forecasts for the neu¬ 
trino mass constraints that will be achievable with mea¬ 
surements of the 21cm power spectrum from future in¬ 
tensity mapping experiments. We consider two surveys, 
both with Phase 1 of the Square Kilometre Array (SKA): 
a wide and deep survey at low redshift (0 < z < 3), us¬ 
ing the SKAl-MID array, and a narrow and deep sur¬ 
vey at higher redshift (3 < z < 6) using SKAl-LOW. 
The forecasts for the latter use 21cm power spectra mea¬ 
sured from our high-resolution hydrodynamic simula¬ 
tions, which extend into the non-linear regime, while the 
former use only linear modes from power spectra derived 
using halofit (Takahashi et al. 2012). 

In all cases we include a number of nuisance param¬ 
eters related to the method used to model the spatial 
distribution of neutral hydrogen, and introduce conser¬ 
vative priors from contemporary experiments to break 
degeneracies between cosmological parameters. We also 
consider more pessimistic choices of instrumental param¬ 
eters than given in the baseline designs of the SKAl sub¬ 
arrays. The result is a conservative set of forecasts for 
the neutrino mass constraints that one can expect to ob¬ 
tain with the SKA. Nevertheless, a number of systematic 
effects could further degrade the constraints; we assess 
their likely importance at the end of this section. 


mological 21cm brightness temperature power spectrum 
and the instrumental noise, Ptot = P 2 icm + Pn, where 
we have assumed the noise to be Gaussian and uncorre¬ 
lated, as well as uncorrelated with the signal. These are 
reasonable assumptions in the absence of more detailed 
instrumental simulations, and if shot noise is negligible 
(this should be the case; see Santos et al. 2015). 

The form of the noise power spectrum depends on 
the type of radio telescope used to observe the 21cm 
emission. The basic division is between interferome¬ 
ters, which coherently cross-correlate the signals from 
pairs of receivers to observe certain Fourier modes on 
the sky (with mode wavelength corresponding to the in¬ 
verse separation of the receivers), and autocorrelation ex¬ 
periments, which construct sky maps pixel-by-pixel from 
the detected (autocorrelation) signals from individual re¬ 
ceivers for many different pointings. The relative merits 
of the two types are discussed in Bull et al. (2015), where 
models for their noise power spectra are also derived. We 
will simply quote them here. The basic expression is 


Pn Tsys^aiea i rj-2^ 

— = r-i’ll IJ I X, 

ttoti21cm " 


(14) 


where Tsys is the system temperature, S'area is the to¬ 
tal survey area, ftot is the total observing time, J^ 2 icm 
is the rest frame 21cm emission frequency (« 1420 
MHz), and r^, = c(l -I- z)^/H(z). The system temper¬ 
ature is approximately the sum of the sky temperature, 
Tsky « 60K (://300MHz)“^-®, and the instrumental tem¬ 
perature, Tinst- There is an effective beam in the radial 
(frequency) direction due to the frequency channel band¬ 
pass, which we model as 


i?ll = exp 


/ (fc||r^(5^)^ 

\ 16 log 2 


where is the channel bandwidth (assumed to be 100 
kHz), and /cy is the Fourier wavenumber in the radial 
direction. The remaining factor of describes the 

sensitivity as a function of angular scale, and is given by 


5.1. Fisher matrix formalism 

To produce our forecasts, we adopt the Fisher matrix 
formalism for 21cm surveys that was described in Bull 
et al. (2015). The Fisher matrix can be derived from 
a Gaussian expansion of the likelihood for a set of 21cm 
brightness temperature fluctuation maps about a fiducial 
cosmological model. For a single redshift bin, the Fisher 
matrix can be written as 

^ If df'k ^p^'9l0gP21cm(fc) 9l0gP21cm(fc) 

- 2 ] "de, 

( 12 ) 

where {0} is a set of cosmological and astrophysical pa¬ 
rameters to be constrained or marginalised over, and we 
have neglected cosmological evolution within the redshift 
bin. The effective survey volume is given by 


l^ff(^) — k),hys 


P 


21cm 




N 


(13) 


with Fphys the comoving volume of the redshift bin. The 
measured power spectrum is a combination of the cos¬ 


FOV/n((i) (interferometer) 

Tthh (autocorrelation). 

For interferometers, the important factors are the instan¬ 
taneous field of view (FOV) and the baseline density dis¬ 
tribution, n{d), where d is the baseline length (related 
to the transverse Fourier wavenumber by k±r = 27rd/A). 
For autocorrelation, we have assumed a Gaussian beam 
of FWHM 6b for each element of an array with Nd dishes, 
and Nh beams per dish. 

The model for the 21cm (signal) power spectrum de¬ 
pends on the range of redshifts and physical scales in 
question. Our high-resolution simulations, described 
above, cover significantly non-linear scales at z < 5.5, 
but cannot be extended into the linear regime or beyond 
z = 3 due to resolution requirements. For z > 3, we use 
the pseudo-RT 1 model as our fiducial HI model. We take 
the 21cm power spectrum measured from the simulations 
at each redshift, and then extrapolate it to larger scales 
using the halofit matter power spectrum, multiplied by a 
(scale-independent) bias that matches the amplitude of 
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the spectra at A: = 0.2/iMpc“^. These scales are suffi¬ 
ciently linear for the halofit power spectrum to be a good 
approximation, and the bias has become almost scale- 
independent by this k value at all relevant redshifts, as 
shown in Fig. 4. 

At redshifts z < 3, we restrict ourselves to approxi¬ 
mately linear scales only, k < 0.2hMpc“^. The power 
spectra are calculated by multiplying the halofit matter 
power spectrum by a bias function and brightness tem¬ 
perature model derived from the halo-based 1 prescrip¬ 
tion for the function Mhi(M, z); see Section 3.5. 

In all cases, we assume that only the monopole of the 
redshift-space 21cm power spectrum can be recovered. 
In principle, one could use the full redshift-space power 
spectrum, which would also allow constraints on the lin¬ 
ear growth rate to be extracted, but we forego this to 
keep our analysis conservative, and because the higher 
order multipoles extracted from the high redshift simu¬ 
lations are noisier than the monopole, making it more 
difficult to numerically differentiate them reliably. 

We calculate the derivatives of logP 2 icm required by 
Eq. (12) numerically, using central finite differences: 
df/dx « [f{x -1- A) — f{x — A)]/2A, where A is the 
finite difference step size and all other parameters are 
held fixed at their fiducial values. While central differ¬ 
ences require twice as many simulated spectra as forward 
differences, they are more accurate for sufficiently small 
step sizes. The step size chosen for each parameter can 
be found from the list of simulations in Table 1. 

5.2. Survey specifications 

Most current and planned 21cm intensity mapping ex¬ 
periments are optimised either for detecting the baryon 
acoustic oscillations at z ~ 1, or studying the epoch of 
reionization at z > 6. While it is not a dedicated IM ex¬ 
periment, Phase 1 of the Square Kilometre Array (SKA) 
will be able to cover the entire redshift range of inter¬ 
est here - from z = 0 toz~6-so we adopt it as the 
reference experiment for our forecasts. 

SKAl will consist of two sub-arrays: SKAl-MID, a 
mid-frequency dish array to be situated in South Africa, 
and SKAl-LOW, an array of low-frequency dipole anten¬ 
nas that will be constructed in Western Australia. SKAl- 
MID will support receivers that cover several different 
bands, spanning the frequency range from 350 MHz up 
to ~ 14 GHz. Only Bands 1 and 2 are relevant for detect¬ 
ing redshift 21cm emission; these are currently planned 
to cover the ranges 350 — 1050 MHz and 900 — 1670 MHz 
respectively. For simplicity, however, we will assume a 
single band from 350 — 1420 MHz with the specifications 
of Band 1 (summarized in Table 4) . 

The MID array is expected to consist of approximately 
200 dishes of diameter 15m, primarily designed for use 
as an interferometer. It will be able to perform IM sur¬ 
veys much more efficiently if used in an autocorrelation 
mode, however (assuming that technical issues such as 
the effects of correlated noise can be mitigated; see Bull 
et al. 2015; Santos et al. 2015), so we will assume an 
autocorrelation configuration in our forecasts. With this 
setup, a total survey area of ~ 25, 000 deg^ should be 
achievable for a 10,000 hour survey; we take these values 
as our defaults in what follows, but also study the effect 
of changing Sarea in Section 5.7. 

SKAl-LOW will cover the band 50 < v < 350 MHz 



SKAl-LOW 

SKAl-MID 

Tinst [K] 

40 + O.lTgky 

28 

Nd X Nt 

911 X 3 

190 X 1 

Umin [MHz] 

210 

375 

t^max [MHz] 

375 

1420 


925 

140 

Sarea [deg^j 

20 

25,000 

hot [hrs] 

10,000 

10,000 


2.75, 3.25, 3.75, 

0, 0.125, 0.375, 0.625, 

.2 bin edges 

4.25, 4.75, 5.25, 

0.875, 1.125, 1.375, 


5.75 

1.625, 1.875, 2.2, 2.8 


TABLE 4 

Array specifications assumed in our forecasts. More 

DETAILED SPECIFICATIONS ARE GIVEN IN SANTOS ET AL. (2015). 


(3 < z < 27), making it one of the few proposed IM 
experiments that overlaps with the redshifts of our sim¬ 
ulations, 3 < z < 6. (The Murchison Wide Field array® 
also partially covers this range, 80 < < 300 MHz, 

3.7 < z < 16.8.) We use only the high frequency half of 
the band, and modify the maximum frequency slightly to 
375 MHz in our forecasts, to better fit with the redshift 
coverage of our simulations. We use a baseline distribu¬ 
tion based on the description in Dewdney et al. (2013). 
All other specifications are given in Table 4. 

The instantaneous field of view of SKAl-LOW is 2.7 
deg^ at z = 3, increasing to 6 deg^ at z = 5. As our de¬ 
fault, we consider a deep 10,000 hour survey over 20 deg^. 
This is comparable to what will be required to detect 
the 21cm power spectrum from the epoch of reionization 
(EoR). A proposed ‘deep’ EoR survey with SKAl-LOW 
will perform 1000 hour integrations in the band 50 — 200 
MHz at five separate pointings on the sky, for exam¬ 
ple, with each pointing covering ~ 20 deg^ at the refer¬ 
ence frequency 100 MHz (Koopmans et al. 2015). These 
pointings cover a total survey area of ~ 20 deg^ at the 
centre of our chosen SKAl-LOW band (« 290 MHz), so 
the surveys could be performed ‘commensally’. We as¬ 
sume double the total survey time, however; while the 
EoR survey may be limited to observing only on winter 
nights to mitigate ionospheric and foreground contami¬ 
nation (Koopmans et al. 2015), these effects are less of a 
restriction at higher frequencies. The dependence of our 
results on S'area and ftot is studied in Section 5.7. 

5.3. Model uncertainties and nuisance parameters 

As mentioned above, the neutrino mass constraints de¬ 
pend on being able to accurately measure the shape and 
amplitude of the 21cm power spectrum, also as a func¬ 
tion of redshift. While the power spectra derived from 
our simulations are self-consistent given a particular HI 
model, most models are calibrated off current, imperfect 
data, and can be inconsistent with one another. Our 
forecasts are therefore subject to a number of modelling 
uncertainties, which we attempt to take into account in 
this section. (We also explore the dependence of our 
forecasts on the assumed fiducial model in Section 5.6.) 

An important but poorly-constrained contribution to 

® http://www.niwatelescope.org/ 
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Massive neutrino datasets 

a{M^) / eV (95% CL) 


-b Planck CMB 

-b Planck CMB 

-b Spectro-z 

Planck Ml/ 

— 

0.461 

0.094 

SKAl-LOW 

0.311 

0.208 

0.118 

SKAl-MID 

0.268 

0.190 

0.104 

SKAl-LOW -b SKAl-MID 

0.183 

0.145 

0.082 

SKAl-LOW -b Planck M„ 

— 

0.089 

0.076 

SKAl-MID -b Planck M„ 

— 

0.071 

0.065 

SKAl-LOW -b SKAl-MID -b Planck M„ 

— 

0.067 

0.058 


TABLE 5 

Marginal 2a (95% CL) constraints on the neutrino mass, for various combinations of surveys and prior information. 


the HI power spectrum is the HI bias, which is a func¬ 
tion of both redshift and scale. To take into account the 
uncertainty associated with the bias model, we incorpo¬ 
rate a simple template-based bias parametrisation into 
our signal model. This is constructed by first measuring 
the bias from the simulations as a function of scale and 
redshift, b(z, k). We then introduce an amplitude, A, and 
shift parameter, a, into the definition of the monopole of 
the redshift-space 21cm power spectrum in each redshift 
bin, such that 


= {zi)[Aib^d{zi,aik)f (15) 

1 + + g/^fid P^{Zi, k). 


We marginalise over the amplitudes in all of our fore¬ 
casts, but the shift parameters are only marginalised for 
the high-redshift survey; the low-redshift forecasts use 
only linear scales by construction, where the bias has no 
scale-dependence, so the {at} are unconstrained. This 
marginalisation procedure is pessimistic in that it as¬ 
sumes no prior constraints on the redshift evolution of the 
HI bias; Ai is left completely free in each bin, while there 
are in fact some existing constraints on its value. The Ai 
parameter also absorbs the uncertainty in other parame¬ 
ters that affect the normalization, such as P, which we do 
not attempt to constrain separately. We have assumed 
somewhat more about the scale dependence - the overall 
shape of the function is fixed - but by allowing the shift 
parameter (and thus the scale at which scale-dependence 
kicks in) to be free in each redshift bin, we are still being 
reasonably conservative with our model. 


5.4. Prior information 

Without additional information about the cosmologi¬ 
cal model parameters, the Fisher matrix can suffer from 
degeneracies, severely degrading the forecast constraints. 
As well as presenting results for the SKAl-LOW and 
MID surveys alone, we also include Fisher matrix pri¬ 
ors for two datasets that will be available contemporary 
with SKAl: a full-sky CMB experiment (Planck), and a 
future spectroscopic galaxy redshift survey (e.g. Euclid 
or DESI). 


Eor the CMB prior, we used an approximate Fisher 
matrix based on the 2015 Planck-only temperature -I- 
polarisation constraints on a two-parameter {Mi, + A"eff) 
extension to ACDM.® To construct it, we took the public 
Planck MCMC chains for this model, calculated the co- 
variance matrix for all parameters of interest, and then 
inverted it. This procedure discards the non-Gaussian 
part of the posterior distribution, which enforces consis¬ 
tency with the Gaussianity assumption of Fisher fore¬ 
casting but loses some information compared with the 
actual Planck constraints. Gonsidering the other caveats 
and simplifications inherent in Fisher forecasting, and 
the approximate Gaussianity of most subspaces of the 
posterior distribution, this is a reasonable approxima¬ 
tion. Note that this method results in a prior that has 
been implicitly marginalised over a large number of sys¬ 
tematic and instrumental effects, such as foregrounds and 
calibration errors. We also marginalised over the effective 
number of relativistic degrees of freedom, A^effj assuming 
that only the GMB gives information on this parameter. 

For the galaxy redshift survey prior, we performed sep¬ 
arate Fisher forecasts for a future spectroscopic survey of 
^ 6 X 10^ galaxies over 15,000 deg^ from 0.65 ^ z < 2.05, 
as was also used in Bull et al. (2015). This is a sim¬ 
ilar specification to the planned Euclid mission, and 
has comparable properties to DESI. Following Amendola 
et al. (2013), we took the bias to be b{z) = a/1 -|- z, 
marginalised as a free parameter in each redshift bin of 
width Az = 0.1. We also assumed that the information 
from the broadband shape of the power spectrum and 
redshift-space distortions can be used, and marginalised 
over a non-linear damping of small scales in redshift- 
space, (Tnl- Gonservatively, we assumed that no infor¬ 
mation on the neutrino mass is obtained from the galaxy 
survey. When combining the galaxy survey Fisher matrix 
with the low-redshift SKAl-MID IM Fisher matrix, we 
also assumed that the surveys are independent, although 
in reality they both probe the same underlying matter 
density field. We do not expect the effect of partially- 
overlapping survey volumes to be large (and note that 

® We used the base_iinu_mnu_plikHM_TT_lowTEB MCMC chains; 
is the effective number of relativistic species. 
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the SKAl-LOW survey volume has no overlap). 

With the priors included, we are therefore forecasting 
for the following parameter set: 

He, Ag, Tis^ All/} base {-^25 Q:i}bias + {b,^NL}gal. 

5.5. Neutrino mass constraints 

The results of our forecasts are summarized in Table 
5 and Figs. 5 and 6. All forecasted errors are quoted at 
the 2a (95%) CL. 

For our approximate Planck 2015 Fisher matrix, the 
constraint on M^, is 0.46 eV, which is consistent with 
the actual Planck 2015 temperature + polarisation up¬ 
per limit of 0.49 eV (Planck Collaboration 2015).^° Our 
forecasts suggest that intensity mapping surveys with 
SKAl-MID and LOW will be able to moderately surpass 
the Planck constraints, yielding a{M^) « 0.3 eV without 
the addition of any prior information to break the strong 
correlations that exist between some of the cosmological 
parameters for the IM surveys. This improves to « 0.2 
eV when LOW and MID are combined. Adding a Planck 
prior on the cosmological parameters only (i.e. ignoring 
information on M^, from the CMB data) improves the 
constraint slightly, reaching ^ 0.15 eV for MID + LOW, 
but it is the inclusion of spectroscopic redshift survey 
data that most strongly breaks the cosmological param¬ 
eter degeneracies; for a Planck CMB -|- spectro-z prior, 
we expect a{Mv) « 0.08 eV for SKAl-LOW -I- MID. 

Similar precision can be gained by combining either 
of the IM surveys with the full Planck data (includ¬ 
ing the neutrino mass constraint), however, without the 
need to add the spectroscopic galaxy survey: we obtain 
cr(M^) = 0.067 eV for SKAl-LOW -k MID Planck. 
This is a result of the IM and CMB constraints hav¬ 
ing different correlation directions for the covariances be¬ 
tween Ml, and certain cosmological parameters - Fig. 6 
shows how the IM and CMB data have somewhat com¬ 
plementary correlation directions in the Mi, — erg plane, 
for example, and Fig. 5 shows correlations with other pa¬ 
rameters. This complementarity is also present in other 
types of large-scale structure survey (Carbone et al. 2011; 
Abazajian et al. 2011). Adding the spectroscopic sur¬ 
vey data improves the constraints on other cosmological 
parameters somewhat, resulting in a slight further im¬ 
provement to a{M^) = 0.058 eV for the combination 
of SKAl-LOW -I- MID -I- Planck. For a best-fit value 
of Ml, — 0.06 eV (i.e. at the minimum bound), this 
final figure would be sufficient for a > 2cr detection of 
non-zero Mi, from cosmological data alone. This would 
be a significant improvement over the 95% upper limit 
of Ml, < 0.17 eV from the Planck 2015 release, which 
combines CMB temperature and polarisation data with 
a compilation of BAO constraints (Planck Collaboration 
2015),^^ and the best current upper limit of Mi, < 0.12 
eV, from Planck CMB + BAO -I- Lya data (Palanque- 
Delabrouille et al. 2015a). Note that an upper limit of 
Ml, < 0.095 eV would also be enough to rule out the 
inverted hierarchy (Lesgourgues & Pastor 2006). 

These figures can be compared with pre-data release forecasts 
for Planck (e.g. Kaplinghat et al. 2003; Perotto et al. 2006; Kitching 
et al. 2008), which predicted u{Mi,) ~ 0.3 — 0.9 eV (95% CL). 

These are the Planck TT,TE,EE+lowP+BA0 constraints; see 
Eq. (54d) of Planck Collaboration (2015). 


It has been shown before that various large-scale struc¬ 
ture surveys, including IM surveys, can provide strong 
constraints on (e.g. Pritchard & Pierpaoli 2008; Mao 
et al. 2008; Carbone et al. 2011; Abazajian et al. 2011; 
Audren et al. 2013; Font-Ribera et al. 2014; Pritchard 
et al. 2015; Sartoris et al. 2015). For example, Pritchard 
& Pierpaoli (2008) found a similar constraint, a{M,^) = 
0.075 eV at Itr, when forecasting for Planck + SKA at 
z > 8 (during the epoch of reionization). Font-Ribera 
et al. (2014) forecast a significantly stronger constraint 
from Planck -|- DESI at lower redshift of a{M,^) = 0.024 
eV, improving to 0.011 eV when Euclid, LSST, and 
Lyman-a forest data are also included. These stud¬ 
ies all used different intrumental and modelling assump¬ 
tions however, and marginalized over different cosmolog¬ 
ical and nuisance parameters, so a direct comparison is 
not possible. To put our result into context, we there¬ 
fore also produced Fisher forecasts for the spectroscopic 
survey with Mi, now included as a parameter, finding 
a{M,y) = 0.060 eV for the combination of spectro-z + 
Planck Mu. Up to the same forecasting assumptions, 
the combination of IM surveys here should therefore be 
competitive with future spectroscopic surveys like Euclid 
and DESI in terms of a neutrino mass measurement. 

Perhaps more significant is the forecast of a{Mu) = 
0.089 eV for SKAl-LOW -I- Planck. While not quite 
enough for a 2a detection of M^ = 0.06 eV, it is neverthe¬ 
less a strong constraint, using a completely different red¬ 
shift range and observation technique to other planned 
surveys. This is despite our pessimistic analysis, which 
marginalises over the amplitude and scale-dependence of 
the bias in each redshift bin, and uses only the monopole 
of the redshift-space power spectrum. The redshift range 
we assumed for LOW (z ^ 3 — 6) has several advan¬ 
tages for matter power spectrum measurements - non¬ 
linear effects are important only on significantly smaller 
scales than at lower redshifts, and radiative transfer pro¬ 
cesses that affect the 21cm signal at higher redshift, in 
the EoR (Semelin & Iliev 2015), are not present. An IM 
survey, piggy-backed on the deep EoR survey that will 
be performed by SKAl-LOW anyway, may therefore be 
an interesting prospect for a more ‘systematic-tolerant’ 
neutrino mass measurement survey (although foreground 
contamination and instrumental systematics are still seri¬ 
ous issues that would need to be resolved; see e.g. Alonso 
et al. 2015; Chapman et al. 2015; Wolz et al. 2015). 

5.6. Dependence on HI prescription and bias 

The various methods used to model the spatial distri¬ 
bution of HI discussed in Section 3 give somewhat dif¬ 
ferent predictions for the amplitude of the 21cm power 
spectrum, especially at high redshift. Clearly this can 
impact the Mi, constraints by (e.g.) changing the SNR 
of the power spectrum detection, so in this section we 
investigate the senstivity of a{Mu) to this choice. We 
concentrate on the SKAl-LOW -I- Planck constraints, 
including the Planck Mi, information. There is a spread 
of a{Mu) values (at 95% CL), ranging from 0.089 eV for 
the default pseudo-RT 1 method, to 0.097 eV (pseudo- 
RT 2), 0.080 eV (halo-based 1), and 0.083 eV (halo-based 
2); corresponding forecasts for M^ vs. as for each HI pre¬ 
scription are shown in Eig. 7. 

At low redshift, uncertainty in the bias amplitude 
parameters is a potential limiting factor for the neu- 
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Fig. 5.— Forecast 2cr marginal constraints on M^, and several cosmological parameters, for various combinations of surveys. SKAl-MID 
constrains the cosmological parameters significantly better than LOW, but the neutrino mass constraints are comparable. 
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Fig. 6 .— Forecast marginal 1- and 2cr constraints on Mj^ and 
erg, for SKAl-LOW and Planck, and for the combination of all ex¬ 
periments considered here. The combined constraint (red) is much 
better than the individual constraints because multiple parame¬ 
ter degeneracies are broken by combining the datasets. The lower 
bound, CMB-hBAOH-Lya, and Planck 2015 95% limits are shown 
as vertical dashed lines from left to right respectively. 


Fig. 7.— Forecast 2a constraints on and as for SKAl-LOW 
-h Planck, for the different methods used to model the HI (see 
Sec. 3). The largest ellipses are obtained by employing the pseudo- 
RT 1 method (green, solid line) and the pseudo-RT 2 method (grey, 
dashed line); the smallest are for the halo-based 1 method (red, 
solid) and halo-based 2 method (blue, dashed). 


trino mass constraint. Without any bias priors, we find 
that SKAl-MID (plus the Planck prior) should achieve 
a{M^) = 0.071 eV. A 1% prior on the bias in all red- 
shift bins improves this slightly to = 0.066 eV, 

while a 0.1% prior further reduces it to 0.058 eV. A high- 
precision bias model is therefore of some use in improving 
the neutrino mass constraints at these redshifts, although 
the gains are too small to justify the difficulty of reaching 
this level of accuracy in practise. At high redshifts, the 
constraints are also relatively insensitive to the bias pri¬ 


ors - for the combination SKAl-LOW -I- Planck, a{Mv) 
improves to 0.082 eV when a 1% prior is applied to the 
Ai parameters. An additional 1% prior on the Ui param¬ 
eters yields a{M^) = 0.079 eV. 

The low redshift forecasts use 21cm power spectra cal¬ 
culated from the Bagla et al. (2010) HI model for the 
function Miii{M,z), which has two free parameters - 
fmin and Umax “ as described in Section 3.5. These ef¬ 
fectively set the normalisation of the 21cm power spec¬ 
trum, and are constrained by the need to reproduce the 
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Fig. 8. — Forecast 2cr constraints on Mu as a function of survey 
area, for the SKAl-LOW and SKAl-MID surveys combined with 
Planck. The current best upper and lower limits on the sum of 
the neutrino masses are shown as dashed lines. The fiducial survey 
areas are marked as crosses. 

observed HI density evolution. They are nevertheless 
subject to some uncertainty, so we also check the effect 
of allowing them to be free parameters. We find that 
there are strong correlations between Vmin and Umax and 
the bias parameters, to the point where a weak prior on 
one of the two (e.g. cr(vmax) ~ 100 kms“^) is needed to 
make the Fisher matrix invertible for SKAl-MID alone. 
Once this has been applied, however, there is essentially 
no correlation between the Bagla parameters and M^, 
meaning that their effect on a{M^) is limited to chang¬ 
ing the signal-to-noise ratio of the 21cm power spectrum 
detection. 

5.7. Dependence on survey parameters 

Finally, we check the dependence of our results on the 
assumed IM survey parameters. Fig. 8 shows a{M^) as a 
function of the total survey area, ^area, for SKAl-LOW 
and MID, combined with the full Planck Fisher matrix. 
Fig. 9 shows the same as a function of survey time, ftot- 

For a fixed survey time of ftot = 10^ hours, the depen¬ 
dence on survey area is relatively mild for both surveys 
- SKAl-LOW yields essentially the same constraints for 
survey areas up to 100 deg^. Similarly, the MID con¬ 
straints improve relatively little above S'area ~ 2000 deg^. 
In both cases the fiducial survey areas are close to opti¬ 
mal for the chosen survey time. 

Fixing the survey areas to their fiducial values, we find 
a greater sensitivity to the chosen ftot- Larger survey 
times than the fiducial value of 10,000 hours are largely 
impractical, but would be required for the constraints for 
either survey to cross the minimum mass value (0.06 eV) 
at the 2 ct level - MID would require 30,000 hours to reach 
this limit, while LOW would require 100,000 hours. MID 
provides a useful improvement over the current best con¬ 
straint for a minimum survey time of ftot ^ 600 hours, 
while LOW requires ftot ^ 3000 hours. Note that all 
survey times should be thought of as ‘effective’ values, 
as foreground cleaning and other effects will increase the 
noise level, requiring an increase in survey time to com¬ 
pensate. 


Fig. 9.— Forecast marginal Icr constraint on Mi, as a function 
of survey time, for the same combination of surveys as in Fig. 8. 
The SKAl-LOW IM survey is taken to be 20 deg^. 

6 . SUMMARY AND CONCLUSIONS 

Neutrinos are one of the most enigmatic particles in 
nature. The standard model of particle physics describe 
them as massless particles, while observations of the so- 
called neutrino oscillations imply that at least two of 
the three neutrino species must be massive. Massive 
neutrinos are therefore one of the clearest indications of 
physics beyond the standard model. As such, one of the 
most important currently-unanswered questions in mod¬ 
ern physics is: what are the masses of the neutrinos? 

Constraints on the neutrino masses arising from lab¬ 
oratory experiments are not very tight, m(7'e)<2.3 eV 
(Kraus et al. 2005), but are expected to shrink down to 
TO(i>e)<0.2 eV in the coming years. On the other hand, 
constraints obtained by using cosmological observables 
are very tight: < 0.12 eV (95% CL) (Palanque- 

Delabrouille et al. 2015a). This constraint is obtained 
by combining data from CMB, BAO, and the Lya-forest 
(Palanque-Delabronille et al. 2013). 

In the near future, 21cm intensity mapping observa¬ 
tions in the post-reionization era will be introduced as a 
powerful new cosmological tool (Bull et al. 2015). These 
observations can also be used to place extremely tight 
constraints on the neutrino masses. In order to extract 
the maximum information from these surveys, one needs 
to understand, from the theory side, the impact that 
massive neutrinos have on the spatial distribution of neu¬ 
tral hydrogen, both at linear and fully non-linear level. 

In this paper we study, for the first time, the detailed 
effects that massive neutrinos have on the neutral hydro¬ 
gen spatial distribution in the post-reionization epoch, 
focusing our attention on their impact on the HI clus¬ 
tering and abundance. We do this by running hydro- 
dynamical simulations with massless and massive neu¬ 
trinos that cover the redshift range 3 ^ ^ ^ 5.5. The 
output of our simulations is corrected a posteriori to ac¬ 
count for HI self shielding and the formation of molec¬ 
ular hydrogen. We use four different methods to model 
those processes. In our fiducial method (pseudo-RT 1), 
the HI self-shielding is corrected by employing the fit¬ 
ting formula of Rahmati et al. (2013a), while we use a 
phenomenological relation where the fraction of molec- 
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ular hydrogen depends on the hydrodynamical pressure, 
based on the Blizt & Rosolowsky and THINGS observa¬ 
tions (Blitz & Rosolowsky 2006; Leroy et al. 2008), to 
model the formation of molecular hydrogen. 

We find that neutral hydrogen is more clustered in 
cosmologies with massive neutrinos, although its abun¬ 
dance, f2Hi(^), is lower. These differences increase with 
the sum of the neutrino masses, and are mainly due to the 
impact that massive neutrinos have on the spatial distri¬ 
bution of matter, which is well known (see Appendix A): 
they suppress the amplitude of the matter power spec¬ 
trum, and thus the abundance of dark matter halos, on 
small scales. The reason for the differences in the spa¬ 
tial distribution of neutral hydrogen between the mass¬ 
less/massive neutrino cases is mainly due to the different 
cosmologies; massive neutrinos barely modify the func¬ 
tion Mhi(M, z). That is, for a fixed dark matter halo 
mass, the HI mass in cosmologies with massive/massless 
neutrinos is very similar, although on average it increases 
slightly with the sum of the neutrino masses. 

We find that the Mhi(M, z) function can be well fit¬ 
ted by Mhi(M, z) = (M/Mo)“ in the redshift and mass 
range covered by our simulations. The value of a de¬ 
creases with redshift and, for a fixed redshift, increases 
with the sum of the neutrino masses. This points out 
another well-known effect of massive neutrinos: in cos¬ 
mologies with massive neutrinos, the spatial distribution 
of neutral hydrogen on small scales can be viewed as the 
spatial distribution of HI in the corresponding massless 
neutrino model at a earlier time. 

In terms of the HI column density distribution func¬ 
tion, we find that our fiducial massless neutrinos model 
reproduces observational measurements in the redshift 
range 3 ^ z ^ 5 very well. In cosmologies with massive 
neutrinos we find a deficit in the abundance of DLAs and 
sub-DLAs, which can be explained by the lower value of 
f^Hi(-z) present in those cosmologies. 

As stated above, the HI is more strongly clustered in 
cosmologies with massive neutrinos. The reason is that 
Mhi(M, z) barely changes, but halos of the same mass 
are more clustered in massive neutrino cosmologies; this 
happens because massive neutrinos suppress the abun¬ 
dance of dark matter halos, and therefore halos of the 
same mass will cluster more strongly in cosmologies with 
massive neutrinos, since they are rarer in those models. 

Even though the HI is more clustered in massive neu¬ 
trino cosmologies, we find that the amplitude of the 21cm 
power spectrum is lower in those models with respect to 
the model with massless neutrinos. The reason for this 
is again that is lower in those models. 

We have forecasted the constraints that the future SKA 
telescope will place on the sum of the neutrino masses 
by means of 21cm intensity mapping observations in the 
post-reionization era. We did this using the Fisher ma¬ 
trix formalism, modeling the spatial distribution of neu¬ 
tral hydrogen in 14 different cosmological models from 
redshift z = 0 to z « 6. At redshifts z > 3 we use hydro- 
dynamic simulations to simulate the spatial distribution 
of neutral hydrogen, while at z < 3 we use a simple an¬ 
alytic model based on halofit. 

Our forecasts are summarised in Table 5. We find that 
with 10,000 hours of observations in a deep and narrow 
survey, SKAl-LOW, employing the interferometer mode, 
will be able to measure the sum of the neutrino masses 


with an error a{My) ~ 0.3 eV (95% CL), lower than 
the one obtained from CMB observations alone (Planck 
Collaboration 2015). A slightly tighter constraint can 
be placed with 10,000 hours of observations with SKAl- 
MID, using single-dish (autocorrelation) mode for a wide 
25,000 deg^ survey. 

Combining the SKA results with CMB data, we fore¬ 
cast that the sum of the neutrino masses can be con¬ 
strained with a much smaller error, as various degen¬ 
eracies with cosmological parameters are broken. For 
instance, by using CMB constraints on cosmological pa¬ 
rameters alone (i.e. neglecting information on My from 
the CMB), we find that SKAl-LOW-l-CMB will be able 
to measure the sum of the neutrino masses with an er¬ 
ror a{My) ~ 0.21 eV (95% CL), while the combination 
SKAl-MID -|- CMB will improve this to a{My) ~ 0.19 
eV (95% CL). On the other hand, if we do include the 
information on My contained in the CMB, we find that 
a{My) ~ 0.09,0.07 for SKAl-LOW -k CMB and SKAl- 
MID -I- CMB, respectively. 

Finally, by also including a prior on the cosmological 
parameters from a spectroscopic redshift survey such as 
Euclid or DESI, we see a modest further improvement. 
For instance, by combining data from the CMB (neglect¬ 
ing information on the sum of the neutrino masses), a 
spectro-z survey, and the SKA (either SKAl-LOW or 
SKAl-MID), we find that a{My) ~ 0.1 eV. As expected, 
the tighter constraint on the sum of the neutrino masses 
can be obtained by combining the full data (including the 
CMB My constraint), a spectro-z survey, and data from 
both SKAl-LOW and SKAl-MID: a{My) = 0.058 (95% 
CL), which represents a significant improvement over the 
current tightest bound arising from CMB-|-BAO-|-Lya- 
forest data of My < 0.12 eV (Palanque-Delabrouille et al. 
2015a). 

In order to check the robustness of our forecasts 
for SKAl-LOW, which were obtained by running high- 
resolution hydrodynamic simulations and modeling the 
HI using the pseudo-RT 1 method, we have repeated the 
analysis by modeling the spatial distribution of neutral 
hydrogen using 3 completely different methods. We find 
that our results are very stable against the model used 
to assign the HI, as can be seen in Fig. 7. 

Our forecasts depend only weakly on the assumed sur¬ 
vey area (see Fig. 8), but are more sensitive to the total 
observation time available (Fig. 9). A 10,000 hour sur¬ 
vey performed ‘commensally’ with a deep EoR survey 
over 5 pointings on the sky with SKAl-LOW is a realis¬ 
tic prospect however, and should be able to put strong 
constraints on the sum of neutrino masses. 

We conclude that massive neutrinos imprint distinc¬ 
tive signatures on the 21cm power spectrum in the post¬ 
reionization era, and that future 21cm intensity mapping 
surveys with (e.g.) the SKA will be able to place very 
tight constraints on the sum of the neutrino masses. 
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APPENDIX 

IMPACT OF NEUTRINOS ON THE TOTAL MATTER DISTRIBUTION 

Here we briefly discuss the impact of massive neutrinos on the total matter spatial distribution (see Fig. 1 for a visual 
comparison of the spatial distribution of total matter, i.e. CDM-|-baryons-|-neutrinos-|-stars, between cosmologies with 
massless/massive neutrinos). We note that this has already been studied in many different works (Ma & Bertschinger 
1994, 1995; Lesgourgues & Pastor 2006; Saito et al. 2008; Brandbyge et al. 2008; Wong 2008; Saito et al. 2009; 
Brandbyge et al. 2010; Rossi et al. 2014; LoVerde 2014; Agarwal & Feldman 2011; Bird et al. 2012; Wagner et al. 
2012; Viel et al. 2010; Villaescusa-Navarro et al. 2013b; Ali-Hai'moud & Bird 2013; Lesgourgues et al. 2013; Bias 
et al. 2014; Massara et al. 2014; Inman et al. 2015; Fiihrer & Wong 2015; Peloso et al. 2015), but remark that the 
majority of N-body simulations with massive neutrinos are not hydrodynamic, with the exception of Viel et al. (2010); 
Villaescusa-Navarro et al. (2013b); Rossi et al. (2014), and one simulation in Bird et al. (2012). 

In Fig. 10 we show the total matter power spectra for the cosmologies with Mi, = 0.0, 0.3 and 0.6 eV neutrinos at 
0 = 3, z = 4 and z = 5, i.e. using the simulations zz+ and We find that massive neutrinos induce a suppression 
in the amplitude of the total matter power spectrum. The suppression increases with the sum of the neutrino masses 
and is almost redshift-independent. The suppression has its physical origin in the fact that the large thermal velocities 
of the neutrinos prevents their clustering on small scales. The bottom panels of that figure display the ratio between 
the matter power spectrum of the massive neutrinos to the massless neutrino model. We obtain the typical shape 
induced by massive neutrinos on that ratio, which can be explained by the extension of the halo model presented in 
Massara et al. (2014). 

In Fig. II we show the halo mass function of the cosmological models with massless and massive neutrinos at 
redshifts z = 3, z = 4, and z = 5. We find that massive neutrinos suppress the abundance of dark matter halos, with 
the suppression increasing with the sum of the neutrino masses, with the halo mass, and with redshift. The weaker 
clustering of matter in cosmologies with massive neutrinos, with respect to their massless neutrino counterpart, is also 
what induces the suppression in the abundance of dark matter halos. 
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Fig. 10.— Impact of massive neutrinos on the total matter power spectrum obtained from our high-resolution hydrodynamic simulations. 
Black lines show the matter power spectrum for the model with massless neutrinos at 2 = 3 (left), z = 4 (middle) and z = 5 (right). 
Magenta and green lines show results for the cosmologies with Miy = 0.3 eV and = 0.6 eV neutrinos, respectively. The bottom panels 
display the ratio between the matter power spectrum of the cosmologies with massive neutrinos to the massless neutrinos model. The 
vertical lines display the Nyquist frequency value of the grid used to measure the power spectrum. 



Fig. 11. — Impact of massive neutrinos on the halo mass function obtained from our high-resolution hydrodynamic simulations. We show 
halo mass function for the model with massless neutrinos (black) and for the models with = 0.3 eV (magenta) and Miy = 0.6 eV (green) 
neutrinos at z = 3 (left), z = 4 (middle) and z = 5 (right). The error bars represent the uncertainty in the halo mass function assuming a 
Poisson distribution. 

















































